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INTRODUCTION 


Human  response  to  impact  environments  is  investigated  using 
instrumented  anthropomorphic  dummies ,  animal  surrogates,  cadavers, 
and  human  volunteers.  This  report  addresses  a  new  technique  for 
analyzing  the  data  from  such  experiments. 

In  a  sled  impact  experiment,  the  dummy,  animal,  cadaver,  or  human 
subject  is  seated  on  a  sled,  which  is  given  a  sudden  impetus  down  a  track 
by  a  piston.  Traditional  observables  used  to  analyze  the  resulting 
motion  include  accelerometer  and  tachometer  velocity  measurements  of 
the  motion  of  the  sled  relative  to  the  track,  photographic  measurements 
of  the  motion  of  fiducials  attached  to  the  subject  relative  to  the  sled, 
and  the  output  of  linear  piezoresistive  accelerometers  attached  to  the 
subject's  body. 

Because  of  the  importance  of  separating  translational  and 
angular  motions,  linear  accelerometer  arrays  are  employed.  However, 
small  percentage  errors  in  measuring  a  large  linear  accelerometer 
output  could  prove  to  be  a  large  percentage  error  in  the  difference 
signal  between  two  accelerometers,  from  which  angular  acceleration  is 
derived.  Thus,  it  is  desirable  to  have  light-weight  instrumentation 
which  can  measure  angular  acceleration  directly. 


Even  without  angular  accelerometers,  improvements  can  be  made 
in  the  traditional  data  analysis  techniques.  The  accelerometer  array 
data  is  currently  analyzed  separately  from  the  photographic  data,  with 
the  latter  being  used,  at  best,  to  provide  an  ad  hoc  correction  to  the 


trajectory  determined  from  the  accelerometer  array  data.  The  more 
modern  technique  is  to  process  all  the  data  simultaneously  through  a 
Kalman  filter  to  obtain  the  trajectory  which  best  fits  the  accelerometer 
array,  photographic,  and  sled  accelerometer  and  tachometer  velocity 
data. 

The  situation  is  analagous  to  navigating  a  cruise  missile  by 
simultaneously  combining  inertial  instrument  and  external  landmark 
sighting  data.  In  the  case  of  a  sled  impact  experiment,  the 
photographic  camera  fixed  to  the  sled  is  looking  at  the  subject, 
rather  than  on  the  vehicle  looking  at  an  external  reference. 

In  order  to  simultaneously  combine  the  various  observables  in  a 
sled  impact  experiment,  this  report  develops  a  quaternion-Kalman  filter 
model  of  the  motions  in  such  an  experiment,  and  analyzes  actual 
experimental  data  witl>  this  technique.  Head  motion  is  particularly 
addressed.  Covariance  analyses  are  performed  to  investigate  the 
possible  accuracy  improvement  that  could  be  obtained  with  an  angular 
accelerometer.  The  Fortran  computer  program  written  to  carry  out  this 
work  is  one  of  the  deliverables  of  this  contract. 

1.1  ANALYSIS  TECHNIQUE 

Formulas  are  derived  for  the  quaternion  representation  of  the 
rotational  motion  of  one  reference  frame  relative  to  another,  and  for 
the  corresponding  angular  velocity  and  acceleration.  Unlike  an  Euler 
angle  representation,  the  quaternion  representation  is  valid  for  arbitrary 
angular  motions,  except  that  a  constraint  is  required  to  restrict  the 
four  quaternion  degrees  of  freedom  to  the  three  degrees  of  freedom 
existing  in  a  rotation. 

A  state  space  dynamical  system  model  is  derived  for  an  impact 
sled  experiment.  The  state  variables  in  units  of  centimeters,  seconds, 
and  radians  are: 
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(a)  Translational  position,  velocity,  and  acceleration  of  the 
sled  frame  relative  to  the  track  (three  states) ; 

(b)  Translational  position,  velocity,  and  acceleration  of  an 
accelerometer  array  frame  attached  to  the  head  relative 
to  the  sled  frame  (nine  states) ;  and 

(c)  Quaternion  rotational  position,  velocity,  and  acceleration 
of  the  accelerometer  array  frame  relative  to  the  sled 
frame  (12  states  with  3  constraints) . 

The  state  dynamical  equations  are  that  the  derivative  of  accleration 
(jerk)  is  equal  to  white  noise. 

Formulas  are  derived  for  photographic,  tachometer,  linear  accelero 
meter,  and  ideal  angular  accelerometer  observables  in  terns  of  the  states 
and  for  the  partial  derivatives  of  the  observables  with  respect  to  the 
states.  Kalman  filter  formulas  are  presented  for  the  propagation  of 
the  state  variables  given  the  observables. 

An  experiment  consnences  with  the  subject  initially  at  rest  with 
a  small  covariance  for  the  uncertainty  in  its  state  initial  conditions. 
The  state  is  propagated  to  an  observable  time  by  the  state  dynamical 
equations.  The  state  covariance  is  also  so  propagated,  with  the 
uncertainty  being  incremented  by  the  white  noise  in  the  dynamics.  This 
plant  noise  is  a  measure  as  to  how  much  the  acceleration  can  change 
between  observable  times,  since  the  dynamical  equations  assume  constant 
acceleration  between  observable  times,  except  for  the  white  noise. 

The  state  estimate  and  its  covariance  are  updated  using  the  observables 
and  the  Kalman  filter  formulas.  The  state  is  then  propagated  to  the 
next  observable  time  and  the  process  repeated,  etc. 

The  proper  weighting  of  past  states  and  present  observables  is 
determined  by  the  relative  size  of  the  propagated  state  covariance  and 
the  covariance  of  the  observables.  Different  observable  types  with 
their  own  units  (g,  inch,  ft/s)  can  be  combined,  since  the  weighting 
covariance  in  the  same  units  creates  dimensionless  combining  terms. 


The  accelerometer  observables  are  one  millisecond  apart,  whereas 
the  photographic  observables  are  two  milliseconds  apart.  There  can  be 
different  numbers  of  observables  at  different  observing  times.  In 
particular,  if  a  sample-and-hold  is  not  employed  in  a  real  time  digital 
data  acquisition  system,  then  the  slew  in  observable  time  for  the 
analog-to-digital  conversion  of  analog  piezoresistive  accelerometer  data 
can  be  taken  into  account.  The  accelerometer  measurements  employed  in 
this  report  were  recorded  on  analog  tape  and  played  back  for  conversion 
to  digital  form,  so  the  accelerometer  data  were  at  the  same  time. 

The  constraints  on  the  quaternion  state  variables  are  taken  into 
account  in  two  ways.  First,  at  each  observation  update  time  three 
additional  perfect  quaternion  constraint  observables  for  the  rotational 
position,  velocity,  and  acceleration  are  added  to  the  real  observables . 
These  perfect  constraint  observables  have  zero  error,  and  a  form  of  the 
Kalman  filter  update  formulas  which  allows  zero  observable  covariance 
is  employed.  Second,  after  state  propagation  and  observable  updates, 
the  quaternion  states  are  normalized  with  the  constraint  conditions. 

It  was  attempted  to  employ  a  Kalman  smoother  in  fitting  the 
trajectory  to  the  observables,  where  the  Kalman  smoother  estimate  is 
the  optimal  combination  of  forward  and  backward  Kalman  filter  estimates. 
The  Kalman  filter  is  propagated  forward  in  time  from  the  state  initial 
conditions.  Then  the  Kalman  filter  is  run  backwards  from  the  last 
state  estimate  of  the  forward  filter  with  completely  uncertain 
covariance  for  this  last  point.  It  turned  out  that  the  backward  filter 
did  not  thereby  yield  good  state  estimates,  because  the  measurements 
at  the  first  and  subsequent  time  points  into  the  past  did  not  give 
complete  observability  into  the  states.  Therefore,  it  is  better  to 
take  presmoothed  observables  and  run  the  forward  Kalman  filter  starting 
with  good  initial  conditions,  rather  than  take  raw  observables  and  try 
to  have  the  Kalman  smoother  get  the  smoothing  effect  of  the  past  and 
future  around  each  point. 


Formulas  are  derived  for  both  lateral  and  fore-aft  (eyes  out) 
impact  sled  tests.  A  quaternion  rotation  is  required  for  the  fore-aft 
orientation  of  the  accelerometer  array  frame  relative  to  the  sled  frame 
as  compared  to  the  orientation  for  the  lateral  case.  The  software 
implementing  the  quaternion-Kalman  filter  model  formulas  allows 
covariance  analysis  of  simulated  data  as  well  as  fitting  to  real  data. 

1 . 2  ANALYSIS  SOFTWARE 

Fortran  IV  software  was  written  implementing  the  formulas  of  this 
report  with  extensive  internal  comments.  The  software  was  written,  tested, 
and  run  on  a  large  Amdahl  (IBM  370-type)  computer  in  such  a  way  that 
it  should  be  straightforward  to  convert  it  to  run  on  a  CD C  6600  computer 
at  Wright-Patterson  Air  Force  Base.  Each  subroutine  commences  with  an 

IMPLICIT  REAL  *8  (A-H,  O-Z) 

statement  to  force  double  precision  16  decimal  place  computations  on  the 
Amdahl  computer.  These  statements  would  be  removed  or  commented  out  on 
the  CDC  6600  and  single  precision  14  decimal  place  computations  employed. 
These  and  other  statements  which  are  expected  to  require  change  on  a 
CDC  6600  are  flagged  by  comment  cards  with  asterisks. 

A  list  of  the  31  subroutines  in  the  sled  impact  data  analysis 
program  is  given  in  Table  1-1.  There  is  a  total  of  3,107  lines  of  code, 
including  comments.  The  memory  required  on  the  Amdahl  computer  is  268K 
bytes.  This  could  be  reduced  to  200K  bytes  a  30K  CDC  6600  60  bit  words 
using  overlay. 

A  typical  300  millisecond  impact  Kalman  filter  run  with  data  every 
millisecond  (nine  linear  accelerometer  array,  six  photographic  fiducial 
coordinate,  and  sled  accelerometer  and  tachometer  observables)  required 

2.3  minutes  of  CPU  time  on  the  Amdahl  computer,  which  would  be  somewhat 
greater  on  a  CDC  6600  computer.  Changes  could  be  made  to  the  software 
to  make  it  slightly  more  efficient. 


Table  1-1.  Impact  sled  test  Fortran  subroutines. 


Subroutine 

Name 

Number 

of 

Lines 

Description 

AAMAIN 

105 

Main  Program 

ANGACL 

126 

Process  Angular  Accelerometer  Observables 

DOT 

24 

Calculate  the  Dot  Product  of  Two  3  Vectors 

DSTATE 

81 

Calculate  State  Transition  Jacobian  Partial 
Derivative  Matrix 

ERRWGT 

78 

Alter  Individual  Observation  Errors  Using  Input 

Error  Weights 

FILTB 

328 

Backward  Kalman  Filter  and  Kalman  Smoother 

F1LTF 

274 

Forward  Kalman  Filter 

INPUT 

547 

Read  Input  Control  Parameters  for  Given  Experiment 

KALCOV 

66 

Perform  Kalman  Filter  State  Covariance  Update 

KALGAN 

46 

Calculate  Kalman  Filter  Gain  Matrix 

LINACL 

196 

Process  Linear  Accelerometer  Observables 

MATMLT 

40 

T 

General  Matrix  Multiplication  AB 

MATMUL 

40 

General  Matrix  Multiplication  AB 

OBSERV 

136 

Process  Observations,  Including  Quaternion 
Constraints 

PHOTO 

111 

Process  Photographic  Fiducial  Observables 

PNOISE 

62 

Determine  the  State  Dynamical  Plant  Noise 

QAUTO 

48 

Quaternion  Automorphism 

QCONJ 

25 

Quaternion  Conjugate 

QNORMX 

57 

Normalize  Quaternion  Rotational  Position,  Velocity, 
Acceleration 

QNPART 

84 

Quaternion  Normalization  Partial  Derivatives 

QPROD 

27 

Quaternion  Product 

QSTRNT 

57 

Apply  Quaternion  Constraint  to  State  Vector  and 

Its  Covariance 

SCLINV 

71 

Scale  a  Symmetric  Matrix  Before  Calling  Matrix 
Inversion 

Table  1-1.  Impact  sled  test  Fortran  subroutines.  (Continued) 


Subroutine 

Name 

Number 

of 

Lines 

Description 

STATE 

95 

Propagate  State  and  Calculate  the  State  Transition 

Matrix 

SYMADD 

27 

Symmetric  Matrix  Addition 

SYMINV 

104 

Symmetric  Matrix  Inversion  In  Place 

SZMUL 

40 

Matrix  Multiplication  SZ 

TRACK 

68 

Process  Motion  Relative  To  Track  Observables 

ZSMUL 

40 

Matrix  Multiplication  ZS 

ZSZMLT 

55 

T 

Matrix  Multiplication  Z  SZ 

ZSZMUL 

55 

Matrix  Multiplication  ZSZT 

Total 

3,113 

Very  extensive  tests  were  performed  on  the  software.  The 
consistency  of  the  state  propagation  equations  and  the  Jacobian 
partial  derivatives  was  checked  by  the  difference  method.  The 
same  technique  was  used  for  checking  observable  and  partial  derivative 
formulas  for  sled  accelerometer  and  tachometer  velocity  observables, 
linear  accelerometer  and  ideal  angular  accelerometer  observables, 
photographic  observables ,  and  perfect  quaternion  constraint  observables. 
Matrix  manipulation  subroutines  were  also  checked.  For  example,  the 
product  of  a  matrix  with  its  calculated  inverse  was  compared  with  the 
identity  matrix. 

Successful  checkout  of  the  Kalman  filter  code  was  done  once 
the  correctness  of  these  various  service  subroutines  was  verified. 

The  ultimate  test  was  that  the  Kalman  filter  tracked  real  data. 

A  plotting  routine  was  written  to  display  the  results  of  data 
analyses.  By  its  nature,  such  a  program  is  installation  specific,  so 
an  equivalent  program  would  have  to  be  specially  written  for  the  CDC 


6600.  Similarly,  installation  specific  software  was  written  to  read 
magnetic  and  punched  paper  tapes  containing  impact  sled  test  data; 
equivalent  software  would  also  be  required  on  the  CDC  6600. 

1.3  FITS  TO  DATA 

Two  lateral  and  two  fore-aft  (eyes  out)  impact  sled  tests  were 
analyzed  with  the  Kalman  filter  software.  There  was  a  dummy  and  a 
human  test  for  each  given  type  of  experiment. 

Lightweight  Endevco  2264-200  piezoresistive  linear  accelerometer 
arrays  were  employed.  In  the  lateral  impact  case  there  was  a  nine-accel¬ 
erometer  array  strapped  to  the  side  of  the  head,  and  in  the  fore-aft 
impact  case,  it  was  strapped  to  the  front  of  the  head.  In  the  fore-aft 
experiment  only  three  of  the  nine  outputs  were  used  to  simulate  the  situa¬ 
tion  of  a  three-accelerometer  array  afixed  to  the  teeth  inside  the  mouth. 
This  latter  arrangement  is  starting  to  be  used  because  the  accelerometer 
array  is  thereby  more  firmly  attached  to  the  skull  than  when  it  is  strapped 
to  the  head  over  the  pliable  skin. 

In  the  lateral  impact  tests,  the  post-fit  data  residuals  show 
that  the  Kalman  filter  is  tracking  the  observables,  which  indicates 
that  it  is  mathematically  working.  Some  of  the  residuals  are  20  to 
30  percent  of  the  size  of  the  10  g  magnitude  impact  event.  Slight 
changes  in  the  input  values  of  the  coordinates  and  orientations  of 
the  fiducials  and  accelerometers  relative  to  each  other  could  cause 
a  dramatic  improvement  in  the  fit  performance  of  the  Kalman  filter. 

There  was  a  three-accelerometer  array  at  the  center  of  the 
dummy’s  head  in  the  dummy  lateral  impact  test.  The  residuals  for 
these  observables,  which  were  not  included  in  the  fit,  were  as  large 
as  22  g.  Better  engineering  understanding  of  the  experimental  setup 
should  help  improve  these  Kalman  filter  predictions. 

The  fore-aft  (eyes  out)  impact  tests  had  a  restricted  amount 
of  data.  Only  the  two  coordinates  in  the  impact  plane  were  recorded 
for  the  head  photographic  fiducials,  besides  having  only  three-accelero¬ 
meter  array  output.  The  Kalman  filter  tracked  the  observables  until 
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the  start  of  the  impact  event,  at  which  point  it  diverged  with  some 
residuals  being  greater  than  10 ^  g.  The  cause  of  this  divergence 
requires  further  investigation. 

The  covariance  simulation  mode  of  the  Kalman  filter  software 
was  utilized  on  a  lateral  impact  test  to  determine  the  improvement  to 
be  obtained  with  a  three-axis  angular  accelerometer.  However,  the  real 
advantage  of  an  angular  accelerometer  is  to  counteract  the  effect  of 
piezoresistive  accelerometer  array  dimensional  and  cross  axis  biases. 

Such  bias  states  would  have  to  be  added  to  the  present  Kalman  filter 
to  obtain  meaningful  results  from  a  covariance  simulation. 

1.4  FURTHER  WORK 

The  sled  impact  Kalman  filter  Fortran  IV  software  should  be 
converted  to  run  on  the  Wright-Patterson  Air  Force  Base  (WPAFB)  CDC  6600 
computer.  The  various  tests  used  in  debugging  the  software  should  be 
performed,  and  the  fits  to  data  reported  herein  repeated.  If  changes  are 
made  to  individual  subroutines  at  The  Charles  Stark  Draper  Laboratory,  Inc. 
(CSDL) ,  they  can  be  incorporated  in  the  CDC  6600  version  of  the  program. 

The  day-to-day  analysis  of  on-going  and  past  experiments  is  best 
done  at  WPAFB.  Refinements  to  the  analysis  technique  are  most  expeditiously 
done  at  CSDL. 

It  is  sound  engineering  practice  to  first  apply  the  Kalman  filter 
tool  to  accelerometer  array  data  in  a  more  controlled  environment  than 
exists  in  a  sled  test.  Namely,  a  special  fixture  is  used  to  impart  known 
linear  and  angular  accelerations  to  an  accelerometer  array  for  scale 
factor  determinations.  The  data  from  such  calibration  tests  should  be 
analyzed  with  the  Kalman  filter  software.  States  would  have  to  be  added 
for  scale  factors,  cross  axis  sensitivities,  misalignments,  etc.  Besides 
checking  the  model  and  software  in  a  controlled  situation,  greater 
understanding  and  more  accurate  scale  factors  and  other  coefficients 
will  be  determined  for  the  accelerometer  array. 


A  more  thorough  understanding  of  the  coordinates  and  orientations 
in  the  various  experiments  should  be  developed,  and  the  data  analyses  of 
this  report  repeated  to  see  if  the  fits  can  be  improved.  In  particular, 
it  should  be  determined  how  well  photographic  fiducial  data  combined  with 
three-accelerometer  mouthpiece  array  data  (rather  than  nine-accelerometer 
array  data)  can  describe  head  motion  in  an  impact  experiment. 

The  software  should  be  made  to  handle  other  experimental  orientations 
besides  lateral  and  fore-aft  impacts.  Nonstationary  initial  conditions 
should  be  tried,  such  as  with  a  sled  traveling  down  the  track  and  hitting 
a  barrier,  e.g.,  a  water  brake. 

The  presentation  of  experimental  results  should  be  improved.  For 
example,  there  should  be  software  to  calculate  angular  velocities  and 
accelerations  at  the  center  of  the  head  and  other  points  beside  the 
center  of  the  accelerometer  array  using  the  quaternion  formulation. 

The  motion  of  other  parts  of  the  body  should  be  addressed .  There 
are  whole  body  models  of  greater  or  lesser  complexity  with  the  various 
discrete  parts  attached  to  each  other.  A  Kalman  filter  with  a  greatly 
expanded  number  of  states  could  be  devised  to  simultaneously  process 
observables  of  the  various  body  parts  using  these  body  models  with  the 
constraint  that  the  body  parts  are  attached  to  each  other.  The  Kalman 
filter  technique  could  also  be  applied  to  accelerometer  and  photographic 
observables  of  the  individual  body  parts  without  the  complexity  and 
difficulty  of  considering  whole  body  models. 

If  the  problem  of  applying  the  Kalman  filter  with  whole  body 
models  is  addressed,  then  other  observables  could  be  employed,  such 
as  body  pressure  readings  on  the  seat  or  belt  restraints.  If  an 
angular  accelerometer  is  employed,  then  the  angular  accelerometer  obser¬ 
vable  formulas  and  code  should  be  changed  to  reflect  the  characteristics 
of  the  actual  rather  them  ideal  instrument 

As  the  sled  impact  Kalman  filter  software  matures,  more  formal 
documentation  should  be  generated.  For  the  present,  the  program  is  to  some 
extent  self-documenting  with  extensive  comments  in  each  subroutine. 
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SECTION  2 


QUATERNION  REPRESENTATION  FOR  ROTATIONS 


2.1  QUATERNION  ALGEBRA 

A  quaternion,  a,  is  a  4-vector: 
a  =  (a3,  a2,  a3,  a4) 


(1) 


where  the  a^  are  real  numbers  (i  =  1,  ...  4).  Addition  and  multiplica 


tion  of  quaternions  are  defined  by: 

a  +  b  =  (a3  +  bj,  a2  +  b2,  a3  +  b3,  a4  +  b4) 
a  •  b  =  (a^bj  —  a2b2  —  a3b3  —  a4b4. 


(2) 


ajb2  +  a2bi  +  a3b4  -  a4b3, 
ajb3  +  a3bj  +  a4b2  -  a2b4, 
ajb4  +  a4bj  +  a2b3  -  a3b2) 


(3) 


Note  that  a-b  is  not  necessarily  equal  to  b-a.  The  product  of  a  real 
number  a  and  a  quaternion  a  is: 


aa 


(4) 


—  (aj ,  -a2 »  —a 3 ,  — a4) 


(5) 


(aalf  aa2,  oa3,  aa4) 

The  conjugate  a*  of  a  quaternion  a  is 
a*  =  (a3 ,  a 2,  a3,  a4) * 

so  that 

(ab)*  =  b*a*,  a**  =  a 

The  absolute  value  squared  of  a  quaternion  a  is  the  positive 
real  number: 

2  ,  0,  0,  0,) 


(6) 


=  aa* 


(a 


+  a?*  +  a-a 


+  a^j 


(7) 


The  multiplicative  inverse  a-1  of  a  is: 


since 


(8) 


a-[fj2  =  f^j2'a  =  (1'  °»  °*  0)  (9) 

the  quaternion  identity.  The  above  properties  are  summarized  by 
saying  that  the  quaternions  are  a  noncommutative  number  field,  a  fact 
originally  discovered  by  Hamilton. 

A  quaternion  number  a  is  a  pure  real  number  if  a*  =  a  and  a  pure 
quaternion  (i.e.,  real  part  zero)  if  a*  =  -a.  For  any  nonzero  quaternion 
q  the  transformation 

a  -*■  qaq-1  (10) 

is  called  an  inner  automorphism  of  the  quaternion  number  field.  It 


preserves  arithmetic  operations  in  the  sense  that 

q(a  +  b)q-1  =  (qaq-1)  +  (qbq-1)  (11) 

q(ab)q-1  =  (qaq-1) (qbq-1)  (12) 

|qaq-1|2  =  |a|2  (13) 

If  q’  =  aq,  where  a  is  a  real  number,  then  q1  and  q  define  the 
same  inner  automorphism.  Therefore,  we  can  impose  the  constraint 

qq*  =  1  (14) 

and  the  inner  automorphism  can  be  written  as 

a  -*  qaq*  (15) 


With  the  constraint  Equation  (14)  there  is  still  a  twofold  ambiguity 
in  the  sense  that  q  and  -q  define  the  same  inner  automorphism. 


2.2 


ROTATIONS  IN  3-SPACE 


Let  v  be  a  quaternion  with  zero  real  part,  v*  =  -v.  It  can  be 
regarded  as  a  vector  in  3-space.  If  q  is  any  quaternion  satisfying 
constraint  Equation  (14),  then  qvq*  also  has  zero  real  part.  Thus, 

v  -*■  qvq*  (16) 

represents  a  linear  transformation  into  3-space.  Conversely,  any 
linear  transformation  (or  rotation)  of  3-space  into  3-space  can  be 
represented  by  a  quaternion  inner  automorphism  Equation  (16)  ^ . 

It  is  a  two-valued  representation,  in  that  both  q  and  -q  represent 
the  same  linear  transformation.  The  quaternion  representation  can  be 
thought  of  as  the  square  root  of  a  linear  orthogonal  transformation  of 
3-space. 

The  quaternion  representation  for  a  rotation  only  requires  one 
constraint  Equation  (14)  on  its  four  components.  The  nine  elements 
in  a  3x3  direction  cosine  matrix  require  six  constraints  to  represent 
a  rotation.  The  three  Euler  angle  parameters  represent  a  rotation  with 
no  constraints.  However,  there  is  the  problem  of  gimbal  lock  (a  partial 
derivative  matrix  singularity  for  certain  values  of  the  angles)  and 
having  dynamical  equations  involving  white  noise  with  sines  and  cosines 
of  Euler  angles  is  less  tractible  than  assuming  white  noise  in  the 
quaternion  components. 

Given  an  initial  choice  for  one  of  the  two  possible  quaternions 
representing  the  initial  orientation  of  the  instrument  frame  [say 
(1,  0,  0,  0)],  the  choice  for  time  increments  after  the  initial  time 
is  defined  by  continuity. 

2.3  CORRESPONDENCE  WITH  MATRICES 

An  explicit  expression  for  automorphism  Equation  (15)  is  by 
Equation  (3) : 
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qaq*  *  [(qj2  +  q22  +  q32  +  q42)*l» 

(qi2  +  q22  -  q32  -  q42)*2  +  2(q2q3  -  qiq4>a3 

+  2 (qiq3  +  q2q4>  a4» 

2 ( q i q4  +  q2q3>a2  +  Cqi 2  -  q22  4  q32  “  q42)a3 
+  2(q3q4  -  qjq2)a4. 

2(q2q4  -  qiq3)a2  +  2(qiq2  +  q3q4)a3 
+  (qi2  -  q22  -  q32  -  q42)a4]  (17) 


By  Equations  (7)  and  (14)  ,  a^  is  left  fixed.  In  particular,  as  proved 
above,  a  3-vector  with  a^  =  0  is  mapped  into  a  3-vector. 


By  Equation  (17)  the  transformation  matrix  corresponding  to  the 
vector  transformation  defined  by  a  quaternion  q  is 


(qi2  +  q22  2(q2q3  -  qjq4)  2(qjq3  +  q2q4) 

-  q32  -  q42) 


2(qiq4  +  q2qs> 


2(q2q4  -  qiq3) 


(qi 2  *  q22  2 (q3q4  -  q3q2) 

+  q32  -  q42) 

2(qjq2  +  q3q4)  (qi2  -  q22 

-  qa2  +  q42) 


(18) 


As  mentioned  above,  q  and  -q  map  into  the  same  linear  transformation. 

Infinitessimal  rotations  by  angles  60j,  602,  603  about  the 
three  coordinate  axes  yield  a  transformation  matrix: 


1 

663 

-602 

-603 

1 

60  j 

602 

-60 1 

1 
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Therefore,  the  quaternion  corresponding  to  this  infinitessimal  transfor¬ 
mation  is: 


q 


662  60j 

2  '  2 


(20) 


or  -q,  where  the  squares  of  the  small  angles:  in  radians  are  ignored. 

Let  0^  be  a  finite  rotation  angle  about  coordinate  direction  i. 

The  rotation  matrix  and  rotation  quaternion  for  each  such  finite  rotation 
are: 


1 

0 

0 

0 

COS  0j 

sin 

0i 

q  = 

(cos 

-|l,  -sin  -|l ,  0,  0) 

0 

-sin  0j 

cos 

ei_ 

(21) 

COS  02 

0 

sin 

02 

0 

1 

0 

«— ► 

q  = 

(cos 

— ,  0,  -sin  1^-,  0) 

-sin  02 

0 

cos 

02 

(22) 

cos  03 

sin  03 

0 

- 

-sin  03 

cos  03 

0 

q  = 

(cos 

0,  0,  -sin  ^3) 

0 

0 

1 

(23) 

or  -q  in  each  case. 

An  Euler  angle  transformation  matrix  can  be  expressed  as  the 
product  of  three  simple  rotation  matrices  such  as  above.  Namely,  let 
(X,  Y,  Z)  be  one  coordinate  system  and  (X’,  Y’,  Z')  another  with  Euler 
angles  given  in  Figure  2-1: 

ft  =  angle  between  the  X  axis  and  the  intersection 

of  the  (X',  Y')  plane  on  the  (X,  Y)  plane  measured 
along  the  (X,  Y)  plane  (0  <  ft  <  360°) 
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Figure  2-1.  Euler  angles. 


I  =  angle  between  the  (X,  Y)  plane  and  the  (X' ,  Y 
(0  <  I  <  180°) 

w  =  angle  between  the  intersection  of  the  (X',  Y' 
plane  and  the  X'  axis  measured  along  the  (X' , 
(0  <  ui  <  360°) 


' )  plane 

>  on  the  (X,  Y) 
Y')  plane 
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t - r* 


is 


All 

3 

cos 

(2 

cos 

u> 

- 

sin 

£2 

sin 

U' 

cos 

1 

Al2 

= 

sin 

£2 

cos 

0) 

cos 

£2 

sin 

GJ 

cos 

I 

A 

sin 

13 

~ 

w 

sin 

I 

a2  1 

= 

-cos 

(2 

sin 

- 

sin 

£2 

cos 

0> 

cos 

I 

A22 

= 

-sin 
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(25) 


It  results  from  a  rotation  about  Z  by  an  angle  (2,  about  the  new  Y  axis 
by  an  angle  I,  and  about  the  new-new  Z  axis  (=  Z')  by  an  angle  id.  Thus, 
the  quaternion  corresponding  to  the  Euler  angle  transformation  matrix 
(25)  is 


(2 


a. 


(cos  0,  0,  -  sin  -) . 


.  (cos  -r,  0,  -  sin  — ,  0) . 


. (cos  y,  0,  0,  -  sin  y) 


n 


,  I  w  £2  .  I  .  u> 

(cos  —  cos  —  cos  —  -  sm  —  sin  —  sin  — , 


ft  .  I  .  w  .(2.1  ui 

cos  -r  sin  —  sin  —  -  sin  —  sm  —  cos  — , 

»  *  4  2  2  2 


-cos  y  sin  y  cos  —  -  sm  —  sm  y  sm  y. 


(2  I  .  w  (2  I  w. 

-cos  y  cos  y  sm  y  -  sm  —  cos  —  cos  j) 


(26) 
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Consider  a  general  direction  cosine  orthogonal  matrix: 


All 

a12 

a13 

A 

=  a21 

a22 

a2  3 

a3  1 

a32 

a33 

The  orthogonality  condition  is  that  the  inverse  is  equal  to  the 

-l  T 
transpose  A  1  =  A  : 


l  A.,  A..  =  6.  . 

,  lk  ]k  i] 


1  if  i  =  j 
0  if  i  M 


Equation  (18)  gives  the  A  corresponding  to  a  quaternion  q.  Given  the 
constraint  Equations  (7)  and  (14) ,  checking  that  matrix  (18)  satisfies 
condition  (28)  is  straightforward. 

Suppose  a  matrix  (27)  is  given  satisfying  the  conditions  (28) . 
Presuming  that  matrix  (27)  can  be  expressed  in  the  form  (18),  it  follows 
that 


2(qi2 

-  q42) 

= 

All 

+  a22 

(29a) 

2(qi2 

-  q22) 

= 

a22 

+  a33 

(29b) 

2(qi2 

-  q32) 

= 

All 

+  a33 

(29c) 

2(q22 

-  q32) 

= 

All 

_  A22 

(29d) 

2  (q22 

-  q42) 

= 

All 

'  A3  3 

(29e) 

2  (q32 

-  q42) 

= 

a22 

-  a33 

(29f ) 

4q2q3 

= 

Al2 

+  a21 

(29g) 

4q2q4 

= 

Al3 

+  a31 

(29h) 

4q3q4 

= 

a23 

+  a32 

(29i) 

4qiq2 

s 

A32 

-  A23 

(29  j ) 

4qiq3 

= 

An 

-  A31 

(29k) 

4qiq4 

= 

a21 

-  a12 

(291) 
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By  Equations  (29e,  n)  and  (29f,  o) 


2  =  (A23  +  A32)  (An  ~  A33) 

2(A12  +  A2i  "  A23  -  A32) 

=  (A2i  -  Ai2)  (An  ~  A33) 

2 (A32  ~  A23  -  A2i  +  A12) 

_  (A13  +  A31) (A22  ~  A33) 

2 (Aj2  +  A21  -  A13  -  A31) 

=  Jbn  -~_Ai2)  (A22  -  A33) 

2(Ai3  -  A31  -  A2i  +  A32)  (32) 

The  expressions  in  Equations  (30),  (31),  and  (32)  with  the  largest 
denominators  should  be  chosen  to  evaluate  q22,  q32,  and  qt*2.  The 

O 

component  q3  can  then  be  calculated  from  any  of  Equations  (29a,  b,  c) . 

A  choice  of  sign  for  qj  determines  the  signs  of  q2,  q3,  and  q 4  by 
Equations  (29j,  k,  1).  There  should  be  no  inconsistencies  in  deriving 
the  q  corresponding  to  an  A  if  orthogonal  conditions  (28)  hold  and  if  no 
reflections  are  involved. 

If  all  the  denominators  in  Equations  (30),  (31),  and  (32)  were 
zero,  A  would  have  to  be  a  diagonal  matrix  with  ±1  on  the  diagonal. 

There  can  be  none  or  two  -Is,  since  otherwise  a  reflection  is  involved 
and  quaternions  do  not  map  into  reflections,  only  rotations.  The  quaternions 
corresponding  to  the  allowed  diagonal  rotation  matrices  are  given  in 
Equations  (21),  (22),  and  (23)  with  6^  =  0®,  180®. 

For  completeness,  the  Euler  angles  corresponding  to  a  rotation 
matrix  A  are  given: 


1 


sin 

I 

=  J  A3l2  + 

A322 

cos 

I 

=  AnA22  " 

A12A21 

sin 

OJ 

=  A3i/sin 

I 

cos 

LJ 

=  A32/sin 

I 

sin 

(2 

=  A21COS  0) 

-  A22sin  u> 

cos 

(2 

=  A^COS  0) 

-  Aj2sin  oi 

0°  <  I  <  180°  (33) 


0°  <  u)  <  360°  (34) 


0°  <  (2  <  360°  (35) 


2.4  ANGULAR  VELOCITY 

T 

Let  r  =  (ri»  r2 »  £3)  be  body  fixed  column  vector  with  the 
transformation  to  another  coordinate  system  being  defined  by 

s  =  Ar  (36) 

where  A  is  a  rotation  matrix  which  is  a  function  of  time.  Then 


0)2 
-WJ 
0 

(37) 

because  of  the  orthogonality  condition  Equation  (28) . 

Consider  the  corresponding  quaternion  expressions  with 
r  =  (0,  r\,  r2.  r3> : 


T  ds 
'  dt 


T  dA 

A 


W3 

-002 


-W3 

0 

wi 


w  x  r 
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J. 


=  qrq* 


*ds  dq  dq* 

q  dtq  =  q  -  r  +  r  £ 


dt 


dt 


Constraint  Equation  (14)  implies 


|3q*  +  qd£*  =  o 

dtH  4dt 


£fsq*  +  2^2  *l\  qdii* 
dt2  dt  dt  dt2 


=  0 


Let 


(38) 

(39) 


(40) 


(41) 


w  =  (0,  wj  u2  w3)  =  2q*^S  (42) 

The  real  part  of  <d  is  zero  by  Equation  (40).  Then  Equation  (39)  can 
be  written 

q*sq  =  “£ur  -  rw  J  (43) 

which  is  the  Poisson  bracket.  Carring  out  the  indicated  operations 
yields 

oor  -  ru>  J  =  (0,  w2r3  -  w3r2,  w3rj  -  wjr3, 

u2r3  -  u3r2)  (44) 

which  is  the  quaternion  equivalent  to  w  *  r.  Thus 

(i)j  -  wj ,  u)2  =  o)2,  oj3  =  w3  (45) 


and  the  quaternion  Poisson  bracket  corresponds  to  vector  cross  product. 


2 . 5  ANGULAR  ACCELERATION 

Differentiation  of  Equation  (37)  gives  the  matrix  form  for 
angular  acceleration  in  the  body  fixed  frame 


T  A1  s 
dt2 


dgj 

dt 


xr+wxwxr 


(46) 


assuming  that  the  vector  r  is  fixed  in  the  body.  Differentiation 
of  Equation  (43)  gives  the  quaternion  equivalent 


,d2s 

* — -c 

dt2 


1 

dw 

dui 

r  — 

1 

+  — 

0) 

,  cor  -  r<u 

2 

dt  r 

dt 

4 

air  -  rw 


(47) 


SECTION  3 


DYNAMICAL  SYSTEM  EQUATIONS 


3.1  COORDINATE  SYSTEMS 

In  Figure  3-1,  a  human  or  dummy  subject  is  seated  on  a  sled 
which  is  constrained  to  move  down  a  track  for  a  lateral  impact  test. 
Figure  3-2  is  a  similar  representation  of  a  fore-aft  (eyes  out)  impact 
test. 

There  is  a  9-accelerometer  array  on  the  side  of  subject's  head 
in  Figure  3-1,  and  a  3-accelerometer  array  inside  the  subject's  mouth 
in  Figure  3-2.  Let  the  coordinate  system  (X,  Y,  z)  be  fixed  in  the 
linear  accelerometer  array  with  origin  at  the  center  of  the  accelero¬ 
meter  array  such  that 

X  *  positive  forward  out  of  the  subject's  head 

Z  =  positive  towards  the  top  of  the  subject's  head 

Y  =  completes  the  right  hand  system  (positive  towards  the 

subject's  left  side) 

Define  a  coordinate  system  (X,  Y,  2)  fixed  in  the  sled  seat  by 

X  *  forward  out  of  the  seat 

Z  =  positive  up 

Y  ®  completes  the  right  hand  system 


Figure  3-2.  Fore-aft  (eyes  out)  impact  sled  test  coordinates 


Nl 


Initially,  the  (X,  Y,  Z)  axes  are  parallel  to  the  (X,  Y,  Z)  axes  for 
an  upright  subject.  The  origin  of  the  sled  coordinates  could  vary 
between  experiments. 

Finally,  the  laboratory  coordinate  system  (X,  Y,  Z)  fixed  in 
the  track  is 


-Y  =  along  the  track  in  the  direction  of  motion 
Z  =  positive  up 

X  =  completes  the  right  hand  system 

For  a  lateral  impact  test,  the  (X,  Y,  Z)  axes  are  parallel  to  the 
(X,  Y,  Z)  axes.  For  a  fore-aft  (eyes  out)  impact  test,  -Y  is  parallel 
to  -X  and  X  is  parallel  to  -Y.  The  origins  of  the  (X,  Y,  Z)  and  (X,  Y,  Z) 
coordinate  systems  are  assumed  to  coincide  at  time  t0. 

At  time  t,  let 


X  =  X 


Y  =  Y  +  yx(t) 


Z  =  Z 


X  =  Y  +  yj (t) 


Y  =  -X 


Z  =  Z 


lateral  impact  test 


fore-aft  (eyes  out) 
impact  test 


where 

yi(to)  =  0 

(49) 

yi  (t)  >  0,  t  >  t0 

with  to  being  the  time  of  initiation  of  the  experiment  when  the  impact 
is  imparted  to  the  sled  in  Figures  3-1  or  3-2. 
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Let  the  transformation  from  accelerometer  array  coordinates 
(X,  Y,  Z)  to  sled  coordinates  (X,  Y,  z)  be 

(0,  X,  Y,  z)  =  (0,  zj (t) ,  z2(t),  z3(t)) 

+  q (t) (0,  X,  Y,  Z)q(t)*  (50) 

where  the  quaternion  q  =  (qj,  q2,  q3/  q4>  defines  a  rotation  and 
z  =  (0,  zj,  z2,  z3)  defines  a  translation.  The  initial  conditions 

might  be 

z(t0)  =  (0,  z10,  z20,  z30)  (51) 

q ( t o )  -  (1.  0,  0,  0)  (52) 

for  an  upright  initially  stationary  subject. 

3.2  EQUATIONS  OF  STATE 

The  state  variables  which  give  the  dynamical  system  representation 
of  the  impact  experiment  are 

yj  =  position  of  the  sled  frame  relative  to  the 

track  frame  (cm) 

y2  =  velocity  of  the  sled  frame  relative  to  the 

track  frame  (cm/s) 

y3  =  acceleration  of  the  sled  frame  relative  to 

the  track  frame  (cm/s2) 

Zj,  z2/  z3  =  translational  position  of  the  center  of  the 

accelerometer  array  frame  relative  to  the 
sled  frame  (cm) 

Z4,  Z5,  zg  =  translational  velocity  of  the  center  of  the 

accelerometer  array  frame  relative  to  the 
sled  frame  (cm/s) 

27,  zg,  zg  =  translational  acceleration  of  the  center  of 

the  accelerometer  array  frame  relative  to  the 
sled  frame  (cm/s2) 


3l*  32*  33*  34 


quaternion  rotational  position  of  the 
accelerometer  array  frame  relative  to 
the  sled  frame  (dimensionless) 


35*  36*  37*  38  =  quaternion  rotational  velocity  of  the 

accelerometer  array  frame  relative  to 
the  sled  frame  (s-1) 

39*  3lo *  3li*  312  =  quaternion  rotational  acceleration  of 

the  accelerometer  array  frame  relative 
to  the  sled  frame  (s-2) 

The  dynamical  equations  satisfied  by  these  24  state  variables  are 


dyx 

dt 

=  y2  +  wx 

dy2 
3 F* 

=  Y3  +  w2 

dy3 

dt 

=  0  +  w3 

dt 

dz . 

_ l 

dt 

“  zj+3  +  w3+j* 

j  =  i» 

•  •  • 

*  6 

=  0  +  w3+j,  j 

GO 

II 

9 

dqj 

dt 

dqj 

dt 

=  V4+Wl2  +  3' 

j  =  i 

9  •  •  • 

,  8 

=  0  +  wl2+j,  j 

=  9,  • 

•  •  9 

12 

(53 

where  w  = 

(wp  •  •  ■  r  W2i|)  is 

zero  mean 

Gausian  white  noise 

in  the 

dynamics . 

The  constraint  Equations  (7) 

and 

( 14)  imply 

1 


(54) 


JI.JU  II  kffsasaa 


4 

l 

i=l 


ql 


i=l 


qiq^ 


i+H 


0 


4 

£  q,q 


i=l 


1  i+8 


-  2  q 

i=l 


2 

i+'-* 


(55) 


(56) 


The  dynamics  in  Equation  (53)  are  simplistic.  Still,  they 
allow  observations  at  different  times  to  be  tied  together,  and  will 
propagate  an  impact  trajectory  which  best  fits  diverse  observation 
types  (e.g.,  accelerometer  and  photographic)  using  Kahman  filler  and 
smoother  state  estimation  formulas.  The  covariance  of  the  integral  of 
the  noise  process  (Wj,  ...  ,  w24)  in  the  Kalman  gain  matrix  has  to  be 
judiciously  chosen  to  match  the  impact  profile  of  the  specific 
experiment  (see  Section  7.4). 

The  dynamics  of  the  z  and  q  states  defining  the  transformation 
Equation  (50)  between  the  sled  and  accelerometer  array  frames  is 
constrained  by  the  human  body  dynamics.  A  finite  element  model  of 
the  human  body  with  more  states  than  used  above  could  have  more 
realistic  dynamic  equations  than  the  simplistic  model  employed  in 
Equations  (53) . 


3.3  INTEGRATION  OF  THE  STATE  EQUATIONS 

Let  <y x (t  ) ,  ...  ,  qi2<t.))  be  the  estimates  of  the  state 

K 

variables  at  time  t^.  For  time  tg,  these  are  just  the  initial  condition  of 
the  experiment,  which  could  be  Equations  (49) ,  (51) ,  and  (52)  plus 
zero  velocities  and  accelerations  for  an  initially  stationary  subject. 
For  time  t^  they  represent  the  propagation  of  the  state  by  the 
dynamical  equations  with  updates  from  the  observations  between  times 
tj  and  t^.  It  is  presumed  that  the  quaternion  states 

(qi(tk),  ...  ,  <112 1\)) 

satisfy  the  constraint  Equations  (54),  (55),  and  (56). 
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The  state  at  time  t  ^  given  the  state  at  time  t^  ami  assuming 
zero  mean  noise  is  by  Equations  (53) 


yi(tk+ilV 

=  yi(V+  y2(tk)(tk+l  ■  V 

+  2  y3(V  (tk+l  "  V  2 

V2(tk+lltk) 

=  y2(V  +  y3(V(Vi  -  V 

y3(tk+Jtk) 

=  y3<y 

Zj(tk+l>V 

=  zj(V  +  zj+3(tk)(tk+i  -  V 

+  2  Zj+6(tk) ltk+l  tk)2, 

j  =  1,  2,  3 

Zj+3(tk+JV 

Zj+3(tk)  +  Zj+6 (tk} (tk+l  “  V 

j  *  1,  2,  3 

2j+6(tkJV 

=  2j+6(W' 

j  *  1,  2,  3 

VVilV 

=  qj(V  +  qj+4{V(tk+l  "  V 

+  2  qj+8(tk)  (tk+l  ‘  V2' 
j  *  1#  ...»  4 

qj+4 (tk+JV  =  qj+4(tk)  +  qj  +  8(V  (tk+l  ~  V 


3  “  If  •••  #4 

V»‘VllV  ■  Ve'V- 

j  «  1,  ...  ,  4  (57) 


where  bars  ate  put  over  the  q  variables  because  they  are  intermediate 
*-o  the  final  variables  which  satisfy  the  constraints  Equations  (54) , 

(55),  and  (56). 
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3.4  QUATERNION  CONSTRAINTS 

The  (qj,  _  ,  q12)  satisfy  the  constraints  except  for  higher 

order  terms  in  (t^+1  -  t^) .  In  order  to  satisfy  the  constraints 
exactly,  let 

yvJv  = 


‘wiy 


’i'VilV2 

1=1 


j  =  1,  ....  4 


1/2  ' 


(58) 


qj+4(tk+l'tk) 


qj+4(tk+lltk) 


qj+8(tk+lltk) 


i^1qi(tk+l^tk)  qi+4  (tk+l!tk)  J 

•  qj(tk+ilV' 

j  =  1  r  #4 

qj+8(tk+litk) 


(59) 


[ 


1 


i=l 


+  qi+4(tk+lltk) 

i=l 

•  Wi'v- 

j  =  1,  ...»  4 


] 


(60) 


Taylor  series  formulas  can  be  used  to  propagate  the  differential 
equations  solutions  rather  than  more  sophisticated  integration  techniques 
because  the  spacing  between  observation  times  t^  and  fck+l  is  of  the  order 
of  a  millisecond,  besides  the  fact  that  the  simplistic  dynamical  model 
employed  has  exact  Taylor  series  solutions. 
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0 
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0  0 


1 

0 


0 

1 


0  0 
0  0 


0 

0 

0 


0 

0 

0 


(66) 


The  integration  of  Equation  (61)  between  times  t  and  t 

JC  Kt  1 


with  white 


noise  w  set  to  zero  yields 


x(tk+ilV  =  *ltk+i'  V  *(V 

where  the  state  transition  matrix  $  is 

V  -  e*p  t^Vi  -  V1 


(67) 


,?0  nIA(tk+r  V]i 

1=0 


(68) 


By  Equation  (63)  this  can  be  written  as 


* ( fck+l '  V 


eXp[A3X3(Vl  *  V] 


exp[A9X9(tk+i  -  V] 


exp  [A 


12X12 


(tk+l  "  V 


(69) 


where  the  matrix  exponentials  would  be  evaluated  by  an  appropriate 
computer  subroutine. 
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rj"r-/ry: 


3zj+3(tk+lltk) 


3zj+6(tk+lltk) 


3zi(V 


3qj(tk+i»V 

3qi(v 


3qj+4(tk+l'tk) 

3qi(tk) 


3qj+8(tk+lltk) 


*3+8 

3qi(V 


0 

i 

= 

j 

1 

i 

= 

j+3 

(tk+l  "  V 

i 

= 

j+6 

j 

= 

1,2,3 

0 

i 

= 

j 

0 

i 

= 

j+3 

1 

i 

= 

j+6 

j 

= 

1,2,3 

1 

i 

= 

j 

(tk+l  ~  V 

i 

s 

j+4 

5(tk+l  ‘  V2 

i 

= 

j+8 

j 

= 

1,2,3, 4 

0 

i 

= 

j 

1 

i 

= 

j+4 

(tk+l  ‘  tk) 

i 

= 

j+8 

j 

= 

1,2, 3, 4 

0 

i 

= 

j 

0 

i 

= 

j+4 

1 

i 

r= 

j+8 

(75) 


(76) 


(77) 


(78) 


(79) 


with  partial  derivatives  not  included  in  Equation  (71)  through  (79)  being 


zero. 
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If  more  complicated  linear  dynamics  are  added  to  Equations  (61) 
than  those  given  by  Equations  (63)  through  (66) ,  then  the  matrix  exponen 
tial  can  be  used.  For  now,  closed  form  Equations  (71)  through  (79)  can 
be  employed. 

4.2  QUATERNION  CONSTRAINTS 

Normalizing  the  results  as  at  the  end  of  Section  3  makes  the 
state  transition  nonlinear,  which  requires  the  use  of  extended  Kalman 
filtering  and  smoothing.  This  requirement  also  arises  from  nonlinear 
formulas  for  the  observables  in  terms  of  the  states. 

For  the  final  normalized  variables 

3qj(tk+JV  "  3qj(tk+lltk)  3qi(tk+l^tk) 


2VV  i=l  3q.(tk+1itk)  "W 


=  1,  ...  ,  32 


Let  6 . .  be  the  Kronecher  delta 
13 


i  j 


and  let 


4 

_ 


Then  constraint  Equations  (58) ,  (59)  ,  and  (60)  imply  for  j  =■  1 


5 1  •  •  •  t 


For  example,  let 


q  =  q  =  (1,0,  0,  0) 

Then 

3q;j/3qi  =  0  (i,  j  =  1,  ...  ,  4) 

which  cannot  be  allowed  in  Equation  (80) .  Since  the  normalization 
correction  will  be  very  small  and  the  Jacobian  partial  derivative 
matrix  is  only  used  in  propagating  the  covariance  matrix,  it  will 
be  assumed  that 

3q 


12  X  12  identity  matrix 


(9 


SECTION  5 
OBSERVABLES 


The  theoretical  formulas  for  accelerometer  array,  ideal  angular 
accelerometer,  photographic,  track  sensor,  and  quaternion  constraint 
observables  are  derived  in  terms  of  the  states 

y 1 >  V2>  Y3>  Z1 »  •••  *  z9 •  5i»  •••  *  512 


5.1  LINEAR  ACCELEROMETER  ARRAY 
Let 


b.  ,  ,  b.  ,  b 
11  l2  13 


position  coordinates  of  linear  pie zoresi stive 
accelerometer  i  relative  to  the  center  of  the 
accelerometer  array  in  the  accelerometer  array 
frame  (X,  Y,  Z) 

position  coordinates  of  accelerometer  i  relative 
to  the  laboratory  frame  (X,  Y,  Z)  fixed  in  the 
track 


c .  ,  c .  ,  c . 
11  12  13 


c 


il' 


direction  cosines  of  the  input  axis  direction 
of  accelerometer  i  relative  to  the  accelerometer 
array  frame  (X,  Y,  Z) 

these  direction  cosines  in  the  laboratory  frame 
(X,  Y,  z)  fixed  in  the  track 


In  quaternion  notation 


b.  * 
i 

(0, 

bil' 

bi2' 

bi3 

s 

bi  * 

(0, 

Bil' 

bi2 ' 

Si3 

(92) 
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By  Equations  (48)  and  (50)  the  quaternion  transformation  between  the 
accelerometer  array  and  laboratory  frames  is 


bi  =  p[  (0/  zjft),  zr(t),  z3(t))  +  qttJb^qU)  *]p* 


+  (0,  0, 

-  yi (t) , 

0) 

(93) 

where  by 

Equation  (23) 

with  0  3 

U 

o 

• 

VO 

o 

• 

P 

=  (1,  0,  0, 

0) 

lateral  impact  test 

(94a) 

P 

o 

% 

o 

'St 

II 

■-*) 

fore-aft  (eyes  out) 
impact  test 

(94b) 

The  acceleration  of  accelerometer  i  relative  to  the  laboratory 
frame  in  the  laboratory  frame  is 

d2b^ 

- -  =  p  {(0,  z7(t),  z8(t),  z9(t)) 

dt2 

+  JL  Iq { t) b . q ( t)  *]  }  p*  +  (0,  0,  -y3(t)  ,  0)  (95) 

dt2  1 


where 

—  [q(t)b.q(t)*]  =  b.q(t)  * 

dt2  1  dt2  1 


ndq  (t)  u  dq  (t)  * 
dt  i  dt 

*  q(t)b 

1  dt2 

with 

dq(t)  _  /dqi(t)  dq2(t)  dq3(t)  dq<,(t)\ 

dt  \  dt  '  dt  " '  dt  '  dt  f 

*  (qs(t)»  q6^)f  q7(t)*  qa^)* 

d2q(t)  m 
dt2 


(q9(t),  qio(t),  qn(t),  qi2<t)) 


(98) 


The  acceleration  sensed  by  linear  accelerometer  i  is  the  vector 
dot  product  of  its  input  axis  direction  and  the  specific  acceleration 
with  both  being  referred  to  the  same  reference  frame,  and  where  specific 
acceleration  is  total  acceleration  minus  the  acceleration  due  to  gravity. 
Either  the  laboratory  frame  or  the  accelerometer  array  frame  could  be 
used.  The  latter  is  preferable,  since  formulas  are  simplified  for  cross 
axis  error  models  and  angular  accelerometers. 

A  stationary  object  in  the  laboratory  frame  (X,  Y,  z)  has  a 
specific  acceleration  of 

g  =  980.3  cm/s2  (99) 

along  the  Z  axis  in  the  vertical  up  direction,  where  g  is  the  combination 
of  the  accleration  due  to  gravity  and  the  centrifugal  acceleration  due 
to  the  earth's  rotation.  Thus,  the  theoretical  value  of  the  specific 
acceleration  of  accelerometer  i  referred  to  the  accelerometer  array 
frame  is  by  Equations  (95)  and  (96) 


ni  =  (°'  nil'  ni2'  ni3) 


=  qp 


*r>* 


d2b. 


dt 


2  +  (0*  0,  g) 


pg 


=  q* [ (0,  z 7,  zq ,  z9)  +  P*(0,  0,  -y3,  g)p]q 


(100) 


where  it  was  chosen  to  keep  the  result  in  terms  of  the  quaternion  state 
variables  instead  of  the  angula 
accelerometer  i  output  is  then 


variables  instead  of  the  angular  velocity.  The  theoretical  value  of 


where  X  is  a  scale  factor  to  change  the  internal  units  of  the  Kalman 
a 

Filter-Smoother  computer  program  (cm/s2)  to  the  accelerometer  observable 
units  and  where  w.  is  a  measurement  bias. 

i 

A  piezoresistive  accelerometer's  output  is  actually  a  voltage 
which  is  input  to  an  analog- to-digital  converter.  It  is  presumed  that 
a  scale  factor  calibration  has  been  performed,  and  that  the  observed 
value  of  observable  has  been  converted  to  g.  Then  the  parameter 
X  in  Equation  (101)  is 

cl 


X  =  1/980.3 

d 


(102) 


For  an  initially  at  rest  experiment,  the  bias  is 


yi  =  Ui0 


Ci3g 


(103) 


where  u . „  is  the  average  value  of  the  stationary  accelerometer  output 
lO 

eind  where 

(0,  cil#  ci2,  ci3)  «  pq(t0)(O,  c^,  ci2,  ci3)q(t0)*P*  (104) 


5.2 


IDEAL  ANGULAR  ACCELEROMETER 


By  Equation  (93),  the  rotation  quaternion  from  the  (X,  Y,  Z) 
reference  frame  in  which  an  singular  accelerometer  is  fixed  to  the 
(X,  Y,  z)  laboratory  reference  frame  is  pq(t) .  By  Equation  (42)  the 
angular  velocity  quaternion  oj  with  components  in  the  (X,  Y,  Z) 
reference  frame  is 


(0,  wj,  o)2#  0)3) 


2q(t) 


*dqj£) 

dt 


(105) 


Thus,  the  angular  velocity  quaternion  in  the  (X,  Y,  z)  laboratory 
reference  frame  is 


(0,  u>\,  0)2#  0)3) 

2p^5E~)q(t,*p* 


(106) 
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The  inertial  angular  acceleration  quaternion  in  the  (X,  Y,  Z)  reference 
frame  is  then 


5  =  (0,  4i »  ?2 '  ?3) 

=  q ( t) *p*^~pq ( t) 

-  2  2S“L*q<t>l  (1071 

If  ci  =  (0,  ,  ci2#  c^3)  is  an  input  axis  direction  cosine 

quaternion  relative  to  the  (X,  Y,  Z)  frame  of  an  ideal  angular  accelero¬ 
meter,  then  the  theoretical  value  of  the  output  of  the  instrument  is 


3 

a.  =  £  E.c. .  radians/s2 

1  j-1  3  ^ 


(108) 


Scale  factor  and  bias  are  ignored  because  an  ideal  instrument  is  assumed. 
A  three-axis  angular  accelerometer  would  make  three  measurements  relative 
to  three  mutually  orthogonal  input  axes  cj#  •  c3  and  give  complete 
visibility  into  the  angular  acceleration. 


5 . 3  PHOTOGRAPHIC 

Triangulation  analysis  of  the  film  from  movie  cameras  attached 
to  the  sled  yield  coordinates  of  fiducials  relative  to  the  sled  frame 
versus  time.  The  fiducials  are  attached  to  the  subject  at  various 
points,  e.g.,  on  the  bridge  of  the  nose  or  beside  an  eye. 

The  sled  frame  coordinates  that  are  given  for  each  fiducial  versus 
time  are  (-X,  Y,  Z) .  Note  that  the  photographic  data  employs  a  left 
hand  coordinate  system,  which  has  to  be  accounted  for  when  processing 
observations  with  the  Kalman  Filter-Smoother  software,  which  uses  the 
right  hand  coordinate  systems  of  Figures  3-1  and  3-2  for  internal 
computations. 
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Let 


zq  =  (0,  zqx»  z02>  2q  3) 


=  initial  sled  frame  coordinates 
of  the  center  of  the  accelero¬ 
meter  array  (cm) 


f0i  =  (Or  foil»  foi2*  f0i3)  =  initial  sled  frame  coordinates 

of  fiducial  i  (cm) 

The  above  quantities  have  to  be  derived  from  experimental  measurements, 
which  are  given  in  the  left  hand  coordinate  system  and  measured  in  inches. 


e^  =  (0,  e^1#  e^2»  e )  =  fixed  position  coordinates  of 

fiducial  i  relative  to  the 
accelerometer  array  frame  (cm) 

It  is  assumed  that  the  fiducials  and  the  accelerometer  array  are  rigidly 
tied  together,  which  is  not  strictly  true,  because 

(1)  An  accelerometer  array  strapped  to  the  head  can  move 
relative  to  the  head  during  the  experiment,  although 
this  is  less  of  a  problem  for  an  array  rigidly  attached 
to  the  teeth; 

(2)  The  fiducials  are  attached  to  the  skin,  which  can 
deform  during  the  experiment,  except  that  certain  areas, 
such  as  the  bridge  of  the  nose,  deform  less  them  others. 

If  the  sled  and  accelerometer  array  frames  are  parallel  at  the  initial 
time ,  then 


foij  -  Z0j'  j 


1,  2,  3 


qo  *  (qoi*  qo2»  qoa»  qo<*> 

represents  an  initial  quaternion  rotation  from  the  accelerometer  array 
frame  to  the  sled  frame  at  the  initial  time,  then 

e.  -  qo*  (foi  “  zO>qo  (110) 
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Given  the  state 


x  “  (y,  z,  q) 

at  time  t,  the  sled  frame  coordinates  of  fiducial  i  at  that  time  are  by 
Equation  (50) 


f . 
1 


(0,  f.,,  f.  ,  f.  ) 
ll  12  13 


(0,  z  ,  z2,  z3) 
+  qe^q* 


(111) 


where 

q  =  (q^  ...  ,  q^)  and  e,  =  (0,  e^,  e^,  ei3) 

Let 

^il'  ^i2'  ^i3  =  theoretical  values  of  fiducial  i  photographic 

observables 

Then 

(112) 

(113) 

(114) 

where  A,  gives  the  conversion  from  centimeters  to  inches: 
b 


Bil  - 

-Af . 
b  U 

ei2  * 

b  i2 

Bi3  = 

A.  f .  , 
b  i3 

A.  =  1/2.54  *  0.3937 

b 


(115) 


5.4  TRACK  SENSORS 

An  accelerometer  mounted  on  the  sled  measures  the  acceleration 
of  the  sled  relative  to  the  track  in  g.  The  theoretical  value  of 
this  observable  is 

*1  *  V3  {U6) 
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A  tachometer  attached  to  a  wheel  running  along  the  track  measures 
the  velocity  Y2  of  the  sled  relative  to  the  track  in  ft/s.  The  theoretical 
value  of  this  observable  is 

>2  =  Xcy2  (117) 

where  A is  the  conversion  from  cm/s  to  ft/s: 

A  =  1/(2. 54  x  12)  =  1/30.48 

c 

=  0.0328084  (118) 


5.5  QUATERNION  CONSTRAINTS 


The  quaternion  constraints  are  introduced  as  errorless  artificial 
observations  6j,  62 ,  63  in  order  to  have  the  Kalman  filter  update  yield 
quaternion  states  which  satisfy  the  constraints.  By  Equations  (54) , 
(55),  and  (56) 

4 

_  E  qf  (119) 

i=l 

4 

=  I  (120) 

1*1 


-Si 

62 


-S3 


4 

E 

i=l 


qiqi+8 


4 

E 

i*l 


(121) 


The  observed  values  of  these  observables  are  zero  with  zero  error. 


SECTION  6 

PARTIAL  DERIVATIVES  OF  OBSERVABLES 


The  partial  derivatives  of  the  observables  with  respect  to  the 
states  are  required  to  calculate  the  Kalman  filter  gain  matrix. 


6.1  LINEAR  ACCELEROMETER  ARRAY 

Let  (xj,  ...  ,  X24)  be  given  by  Equation  (62).  Differentiation 
of  Equation  (101)  yields 


ax.  3  an.. 

^  =  j!i 


where  the  quaternion  partial  derivatives  are 

/  3nn  9\2  ani3 

a*£  "  V  3x* '  3x* '  3x* 

By  Equation  (100) 


827 


0,  l  -  1|  2 


q*p* (0,  0,  -1,  Q)pq 

0/  &  *  I#  ...  #6 


q* (0,  1,  0,  0) q 


(122) 


(123) 


(124) 


(125) 


(126) 


(127) 
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where,  for  example,  Sq/Sq^  is  a  quaternion  with  zero  entries  except 
for  a  1  at  position  1. 


6.2 


IDEAL  ANGULAR  ACCELEROMETER 


”7  , 


Differentiation  of  Equation  (108)  yields 


3a. 

< 

3 

=  E 

j=3 

~ — UC  .  . 

.  3X£  13 

(133) 

where  the  quaternion  partial  derivatives  are 

iL 

3xa 

_  /n  35i  352  3Ca  \ 

\  '  3x£  3x£  3x£  J 

(134) 

By  Equation 

(107) 

2L 

3y£ 

=  0, 

A  =  1,  2,  3 

(135) 

35 

3*A 

*  0, 

A  =  1,  ...  ,  9 

(136) 

H 
3q) A 

=  2 

3q*  d2q 

3q£  dt 

+ 

3q»  dg  dq*  j*3L 

3q£  dt  dtq  q  dt  dt  3q£ 

£  =  1#  f  4 

(137) 

3£ 

3qA 

=  2 

q*J_ 

q  3q£  \  dt  /  dt  H 

+ 

A  =  5 ,  ...  ,  8 

(138) 

35 

3qi 

.  2q.»  (£n) 

3V  /• 

A  =  9,  ...  .  12  (139) 
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6.3 


PHOTOGRAPHIC 


Differentiation  of  Equations  (112) ,  (113) »  and  (114)  yields 


!!u 

3f£l 

3x* 

"Ab  3x£ 

(140) 

aai2 

x  3fi2 

(141) 

3x* 

b  3x£ 

3x£ 

= 

3f .  , 

X  — ^ 
b  3x£ 

(142) 

where  the  quaternion  partial  derivatives  are 

3f . 

_ l 

3x£ 

= 

/  3fil  3fi2  >fi3) 

V '  3*e '  3-t '  3*« 1 

(143) 

By  Equation 

(111) 

3f . 
* 

s 

0,  l  »  1,  2,  3 

(144) 

3zj 

= 

(0,  1,  0,  0) 

(145) 

3f . 

_ 

3z2 

= 

(0,  0,  1,  0) 

(146) 

3f± 

3z3 

= 

(0,  0,  0,  1) 

(147) 

H  =  4,  ...  »  9 


3q*  1 


+  qe. 


3q* 


If 


l  ”  5 ,  ...»  12 


(148) 


(149) 

(150) 
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6.4  TRACK  SENSORS 

By  Equations  (116)  and  (117) 

3Y 

1 

3y  =  Xb 

3 

ay 

2 

T  *  \ 

3y  c 

2 

with  all  other  partial  derivatives  being  zero. 


(151) 


(152) 


6.5 


QUATERNION  CONSTRAINTS 


Differentiation  of  Equations  (119) ,  (120) ,  and  (121)  yields 


36. 

_ l 

3y„ 


=  o. 


36  . 

_ l 

3z„ 


=  0, 


3  6  j 

3% 


-2q 


£' 


36  j 


=  o. 


362 

"^7 


=  q 


It  4 


36  2 

777 


=  V"' 


362 

aq» 


=  0, 


i  =  1,  2,  3;  £  =  1,  2,  3  (153) 

i  =  1,  2,  3;  t  •  1,  ...  ,  9  (154) 

£  =  1,  ...  ,  4  (155) 

£  =  5,  ...  ,  12  (156) 

£  =  1,  ...  ,  4  (157) 

£  =  5,  ...  ,  8  (158) 

£  =  9,  ...  ,  12  (159) 
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363 

=  q*+e' 


363 


36j 

3qa  =  q£+8' 


i.  =  1,  ....  4 


J l  =  5 ,  ...»  8 


£  =  9,  ...  ,  12 


(160) 


(161) 


(162) 


6.6  CHECK  OF  PARTIAL  DERIVATIVES 


The  coding  of  the  observable  partial  derivatives  can  be  checked 

by  the  difference  method.  Namely,  given  a  state  vector  j 

— ►  j 

X  ~  t  •  .  .  t  X2 4 )  /  ] 

< 

— i y  * 

the  theoretical  value  of  observable  s. (x)  is  computed  and  also  3s. /3x . . 

1  i  3 

The  state  vector  is  changed  by  the  increment  1 


->• 

A. 

1 


(0/  f  0,  ^ 9  •••  *  ^ 


(163) 


xt43  *  lx . xj-r  xj  *  V  xj*i- 


. .  ,  X24) 


(164) 


Then  it  is  checked  that 


3si(x  +  A^) 


3x . 
1 


si  (x) 

3x  . 

1  J 


-+■->■  -+ 
s  (x  +  A  . )  -  s .  (x) 
i  li 

A  . 

1 


(165) 


The  coding  of  the  partial  derivatives  of  normalized  quaternion 
components  with  respect  to  unnormalized  quaternion  components  and  the 
partial  derivatives  of  the  state  components  at  time  t  ^  with  respect 
to  those  at  time  t.  can  also  be  checked  by  this  method. 


SECTION  7 

KALMAN  FILTER  AND  SMOOTHER 


Let  the  state  vector 


(xj  t  • • •  X24 ) 


(166) 


be  defined  by  Equation  (62) ,  where  the  superscript  T  denotes  transpose 
changing  the  row  vector  in  a  column  vector.  Let  the  state  transition 
from  time  t^  to  tfc  be 


*(tk+i*V  =  +  "(tW  V 


k  *  0,  1,  ...  ,  N 


(167) 


where  W  is  the  integral  of  the  zero  mean  Gaussian  white  noise  process 
wf  so  that 


E (W)  =0 

where  E  denotes  expectation  (see  Section  7.4) .  Let 


(168) 


s(tk) 


Q(tk+1»  t^)  =  covariance  matrix  of  W(tk+J,  so  that 

Qi.  =  E (VMT),  i,  j  =  1,  ...  ,  24  (169) 

Let  the  observable  vector  be 

f 

(uj*  •••  1  ®  f  ®i»  •••  >  y  '  y  *  ®i*  62 »  ^3) 

a  D  1  2 

(170) 
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where  its  length  can  vary  at  different  observable  tiroes.  The  observables 
in  terms  of  the  states  are 

*(  W  =  «{*(tk+l))  +  ?k+l  (171) 

where  rk+1  is  zero  mean  Gaussian  noise  with  covariance  matrix  R^+  , 
uncorrelated  from  one  observing  time  to  the  next.  It  can  be  taken  as 
a  diagonal  matrix  with  elements  that  represent  the  squares  of  standard 
deviations  of  the  observable  measurements.  These  standard  deviations 
are  zero  for  the  quaternion  constraint  observables  6j,  62*  63. 

Formulas  are  now  given  for  forward  and  backward  extended  Kalman 

filters  and  for  the  Kalman  smoother  which  is  the  optimal  combination 
121 

of  the  two  filters  . 


7 . 1  FORWARD  FILTER 


At  the  initial  time  t0,  the  state  initial  condition  for  the 
forward  filter  xf(t0)  is  specified  by,  e.g.,  Equations  (49),  (51), 
and  (52)  with  velocities  and  accelerations  zero  for  an  experiment  with 
the  sled  initially  at  rest.  The  covariance  Pf(t0)  of  the  initial  state 
also  has  to  be  specified.  It  can  be  taken  to  be  a  diagonal  matrix  with 
elements  that  represent  the  squares  of  the  standard  deviations  of  the 
uncertainty  of  the  initial  conditions. 

The  expected  value  of  the  state  x^  (t^+1 1 t^)  at  time  t^+1  given 
the  state  x^(t^)  at  time  t^  is  calculated  from  Equation  (167)  with 
W  =  0  in  the  form  of  Equations  (57)  through  (60) .  The  state  transition 
Jacobian  matrix 


9*f(tkJy 

3VV 


3xfi(th+ity 

3VV 


(172) 


is  calculated  from  Equations  (71)  through  (80)  and  (91) .  For  the  first 
12  linear  state  variables  it  is  the  matrix  exponential  state  transition 
matrix.  For  the  last  12  state  variables  the  matrix  exponential  is  modified 
by  the  partial  derivatives  of  the  nonlinear  quaternion  constraints,  which, 
however,  are  assumed  to  form  the  identity  matrix. 
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The  state  covariance  matrix  Pf(tk+i|tk)  at  time  t^+J  given  the 
covariance  P^tt^)  at  time  ^  and  the  covariance  Q(tk+1»  t^  of  the 
dynamical  noise  is^ 

WJV  '  WflVJLi  *  «“*+.■  V  <173> 

T 

where  j  is  the  transpose  of  J. 

The  theoretical  values  of  the  observables  H  are  calculated  from 
Equations  (100)  through  (121) .  The  observables  Jacobian  matrix 

Lk+1  -  —  - - —  (174) 

3xf(tk+ilv 

is  calculated  from  Equations  (122)  through  (162). 

(21 

The  Kalman  filter  gain  matrix  is 


Wt-JV^+l  [Lk+lP£lt3(+lltk)Lk+l 


Note  that  this  form  of  the  gain  matrix  allows  zero  rows  and  columns  in 
R,  as  is  the  case  when  the  perfect  quaternion  constraint  observables 
are  included.  However,  the  total  matrix  in  the  brackets  [ ]  has  to  be 
inver table. 

It  is  possible  that  for  multiple  observations  at  a  given  time 

some  combinations  of  the  observables  would  be  essentially  equivalent, 

T 

which  might  cause  noninvertability  of  the  LP^L  portion  of  the  matrix 
in  the  brackets  [] .  However,  the  nonzero  portion  of  the  observation 
error  matrix  R  will  likely  make  the  sum  in  brackets  []  invertable.  If 
not,  observations  at  a  given  time  can  be  passed  through  the  Kalman 
filter  sequentially  along  with  the  3  quaternion  constraint  observables . 
In  this  manner,  complete  processing  of  n  observations  at  a  given  time 
would  involve  n  inversions  of  4x4  matrices,  instead  of  one  inversion  of 
an  (n+3)x(n+3)  matrix. 


If  the  data  acquisition  system  makes  the  A/D  conversion  measure¬ 
ments  at  distinct  times,  rather  than  using  a  sample  and  hold  to  obtain 
everything  at  the  same  time,  then  this  sequential  processing  would  have 
to  be  done  anyway. 

The  expected  value  of  the  state  x  (t  )  at  time  t  given  the 

t  K+ 1  -  - ,  K+ 1 

state  at  time  t,  and  the  observables  at  time  t.  ,  is 1 
k  k+1 

'  VVilV  +Kk+1 

-  “‘VhJV11  11761 

The  covariance  Pf(t^  +  1)  of  the  state  at  time  t^  given  the  state  at 

time  t,  and  the  observables  at  time  t,  .  is^ 
k  k+1 

VW  '  [I-*WWVVilV 


where  I  is  the  identity  matrix. 

The  above  formulas  define  the  extended  Kalman  forward  filter 
with  the  nominal  trajectory  for  the  extended  filter  linearization  being 
updated  each  observation  point.  Extended  rather  than  linear  Kalman 
filtering  is  required  because  of  the  nonlinearities  in  the  quaternion 
constraints  and  in  the  formulas  for  the  observables  in  terms  of  the 
states. 

The  quaternion  constraint  information  appears  implicitly  in  the 
Kalman  update  Equation  (176) ,  because  the  constraints  are  taken  to  be 
perfect  observations.  However,  the  linearization  employed  in  extended 
Kalman  filtering  could  yield  an  update  (xjttj^),  —  ,  x24<tk+1^  for 
which  constraints  are  not  satisfied  by  the  quaternion  part  (xisO^j), 
...  ,  X24(tjt+1)).  Therefore,  the  normalization  Equations  (58),  (59), 
and  (60)  should  be  applied.  Let  x,  P  be  the  state  and  its  covariance 
matrix  before  the  quaternion  normalization  and  let  x,  P  these  quantities 
after  the  quaternion  normalization  at  time  t^+1 .  Then 


(178) 


3x 

-> 

3x 


12x12 


12x12 


12x12 

3q 

3q 


(179) 


with  the  12x12  matrix  (3q/3q)  being  given  by  the  Equation  (91)  identity 
matrix,  for  the  reasons  given  at  the  end  of  Section  4.2. 


7 . 2  BACKWARD  FILTER 

At  the  last  time  point  t  in  the  data,  let 

n 

Vv  ■  vv  >i8°' 

The  covariance  of  this  initial  value  for  the  backward  Kalman  filter  is 

pK(t  )  =  »,  P.  (t  )_1  =  0  (181) 

d  n  on 

since  the  Kalman  smoother  state  is  the  optimal  combination  of  the  forward 

and  backward  Kalman  filter  states  weighted  by  P,1  and  P,  )  and  the  Kalman 

f  b 

smoother  state  has  to  equal  the  forward  Kalman  filter  state  at  time  tn, 

where  both  contain  all  the  observable  information  from  time  to  to  time 

t  . 
n 

The  expected  value  of  the  state  x^(t^_1|tJc)  at  time  t^  given  the 
state  x^tt^)  at  time  t^  is  calculated  from  Equations  (57)  through  (60) 
replaced  by  t^ 

Jacobian  matrix 


so  that  (t^  -  t^)  is  negative.  The  state 


Jk_,  .  avvJv 

3VV 

Is  calculated  from  Equations  (71)  through  (80)  and  (91) . 


(182) 
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The  state  covariance  matrix  P^(tt_1 | t^)  at  time  t^  given  the 


covariance  P^Ct^)  at  time  t^  and  the  covariance  Q(t^_1 ,  t^)  of  the 
dynamical  noise  (same  in  the  backward  and  forward  directions)  is  ^ 


WA>  ■  (183> 

For  k  =  N,  Pfa(tN)  =  “>,  which  implies  that  P^tt^  It^)  =  ®. 


The  theoretical  values  of  the  observables  H  is  calculated  from 
Equations  (100)  through  (121) .  The  observable  Jacobian  matrix 


’“•Vvi  i  tk>> 


k-X 


(184) 


is  calculated  from  Equations  (122)  through  (162) . 


.  12] 

The  Kalman  filter  gain  matrix  is 


Vi 
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(185) 


The  expected  value  of  the  state  x  (t  )  at  time  t  given  the  state 

D  K~±  r  Q1 

at  time  t^  and  the  observables  at  time  t^  ^  is 1 


-  WjV  +Vi‘*(Vi' 
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The  covariance  P^t^  ^)  of  the  state  at  time  given  the  state  at 

time  t^  and  the  observables  at  time  t,_  .  is^2^ 


Pb(tk-1) 


Vl 
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(187) 


These  forms  of  the  Kalman  filter  equations  are  ideal  with  the  perfect 
quaternion  constraint  observables  for  all  except  the  first  step  in 'the 
backward  filter  process. 

To  go  from  time  t  to  t.,  ,  with  P,  (t  )  =  the  following 
N  N-l  b  N 

alternative  form  of  the  Kalman  filter  formulas  can  be  used1  J : 

Wi>'‘  ■  11881 

Vi  ■  WX-A-,  11891 

i-l  t  -1 

where  P^tt^  1 | tN>  =  0.  However,  the  matrix  j  can  be 

singular  because  certain  states  might  not  be  directly  observable 

from  the  given  measurements  at  time  t  .  Thus ,  P,( t  , | t  )  1 

N-l  b  N-l  N 

can  be  taken  as  a  diagonal  matrix  with  very  small  diagonal  entries. 

The  covariance  for  the  quaternion  constraint  observables  should  also 
be  small  but  nonzero  so  that  R~*  exists. 

It  turns  out  that  for  the  impact  experiments  investigated  in 

this  report,  the  backward  filter  step  from  time  t  to  t  could  not 

N  N-l  _j 

meaningfully  be  made,  neither  with  small  nonzero  entries  in  P,  (t  It.,) 
j  b  N-l '  N 

and  RN_1  Equations  (188)  and  (189) ,  nor  with  large  but  less  than 

infinite  entries  in  P  (t)  rn  Equations  (185)  and  (187).  The  forward 

b  N 

filter  worked  well  because  it  started  from  known  initial  conditions  at 

time  t0.  In  retrospect  it  is  clear  that  the  backward  filter  cannot 

be  employed  starting  with  complete  uncertainty  at  time  t  unless  the 

N 

measurements  at  time  t  give  complete  observability  into  the  states. 

In  this  situation  it  is  better  to  smooth  the  observables  before  passing 
them  through  the  Kalman  forward  filter  rather  than  using  an  optimal 
combination  of  the  forward  and  backward  Kalman  filters  to  smooth  the 
states. 

If  the  backward  filter  were  employed,  the  quaternion  normalization 
described  at  the  end  of  Section  7.1  should  be  applied  to  x^ft^  j)  and 
then  the  smoother  formulas  of  Section  7.3  utilized. 


7.3 


SMOOTHER 


If  both  the  forward  and  backward  Kalman  filters  are  run,  the 

Kalman  smoother  state  estimate  is  the  following  optimal  combination 

[2] 

of  the  two  filter  estimates  : 

P_1  =  p”1  +  P”1  (190) 

f  D 

x  =  Pfp"1  xf  +  p"1  (191) 

where  the  equations  are  evaluated  for  any  time  t  (tQ  <  t  <  t  )  with 
nonsubscripted  quantities  being  the  smoother  results. 

Because  the  quaternion  constraints  are  nonlinear,  x  no  longer 
satisfies  those  constraints  exactly  even  though  xf  and  x^  do.  Therefore, 
the  normalization  of  x  and  the  updating  of  P  must  be  performed  as 
described  at  the  end  of  Section  7.1. 

7.4  PLANT  NOISE 

Successful  application  of  the  algorithms  described  above  requires 
a  choise  of  Q  which  allows  for  the  unknown  but  bounded  variability  in 
jerk  (=  change  in  acceleration) ,  yet  optimally  combines  the  information 
contained  in  the  observables  at  a  given  time  with  the  information 
Contained  in  past  states  (and  also  future  states  if  a  smoother  could 
be  employed) . 

Suppose  in  a  given  impact  experiment  the  observations  are 
uniformly  one  millisecond  apart.  Suppose  the  acceleration  level  of 
the  order  of  lOg  could  change  by  up  to  0.5g  in  a  millisecond.  Then 
the  standard  deviations  of  the  variability  in  the  state  in  one  millisecond 
aure 

acceleration:  500  cm/s2 
velocity:  0.5  cm/s 

_  3 

position:  0.5  x  10  cm 
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The  squares  of  these  values  can  be  taken  as  the  covariances  in  the 
diagonal  Q(t^>  t^+1)  matrix,  except  that  the  small  uncertainty  in  the 
position  covariance  could  be  increased  somewhat  to  account  for  other 
sources  of  uncertainty  besides  integrated  acceleration  uncertainty. 
However,  the  propogated  state  uncertainty  will  likely  give  appropriate 
magnitudes  for  the  uncertainties  in  position  and  velocity  for  closely 
spaced  observations. 


For  the  quaternion  rotation  uncertainties,  appropriate  standard 
deviations  for  changes  in  one  millisecond  could  be 

angular  acceleration:  10  rad/s2 

_2 

angular  velocity:  10  rad/s 

angular  position:  10-5  rad 

for  the  impact  experiments  of  concern  in  this  report. 


7.5  MATRIX  MANIPULATION 

Symmetric  covariance  matrices  should  be  calculated  and  stored 
in  the  computer  in  lower  diagonal  form.  That  is,  a  syrametrix  matrix 
S  is  stored  in  a  vector  array  as 

!S(J  +  (K(K-l) ) /2) ,  J  <  K 

(192) 

S (K  +  (J(J-l) )/2) ,  K  <  J 

The  usual  Fortran  storage  assignment  for  a  general  NXN  matrix  M  is 

M(J,  K)  =  M( J  +  (K-l)N)  (193) 

Subroutines  are  required  to  multiply  a  symmetric  matrix  by  a  nonsymmetric 

T 

matrix  (MS  or  SM)  to  get  a  nonsymmetric  results,  to  compute  MSM  to  get 
a  symmetric  result,  and  to  multiply  general  matrices. 

The  computer  subroutine  to  calculate  the  inverse  of  a  symmetric 
matrix  stored  in  lower  diagonal  form  by  the  Gauss-Jordan  direct  method 
should  choose  ail  the  pivot  elements  on  the  diagonal,  so  no  interchange 
of  rows  or  columns  is  necessary. 
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There  can  be  numerical  problems  with  matrix  entries  being 

badly  out  of  scale  relative  to  one  another,  especially  with  matrix 

inversion.  Clever  choice  of  units  could  help  mitigate  the  problem. 

However,  it  is  easier  to  choose  convenient  units  (cm,  s,  g,  etc.) 

and  then  use  automatic  scaling  before  performing  a  matrix  inversion. 

Namely,  let  S  *  (S  .. )  be  a  symmetric  matrix.  Define  a  diagonal  matrix 
3* 

D  =  (D..  )  with  diagonal  entries 


33 


(194) 


where  the  diagonal  elements  of  S  are  no  doubt  positive  because  it  is 
probably  a  covariance  or  related  matrix.  Define  a  new  symmetric  matrix 

U  =  DSD  (195) 


The  diagonal  elements  of  U  are  1,  and  those  off  the  diagonal  are  likely 
not  way  out  of  scale,  so  that  it  is  numerically  easier  to  invert.  Then 
the  inverse  of  S  is 


1 . 


SECTION  8 
DATA  ANALYSIS 


Two  lateral  and  two  fore-aft  (eyes  out)  impact  sled  tests  were 
analyzed  with  the  software  written  using  the  formulas  presented  in  this 
report.  One  experiment  of  a  given  type  involved  a  dummy  subject  and 
the  other  a  living  human  subject.  A  nine-accelerometer  array  (9-TAP) 
using  Endevco  2264-200  piezoresistive  accelerometers  was  strapped  to 
the  right  side  of  the  subject's  head  for  the  lateral  impact  tests,  and 
a  similar  three-accelerometer  array  was  used  for  the  fore-aft  impact 
tests.  The  piezoresistive  accelerometer  voltage  outputs  were  converted 
to  g  units  using  precalibrated  scale  factors.  Biases  were  determined 
from  the  first  few  data  points  in  each  experiment  when  the  subject  was 
stationary  before  the  impact  event. 

Head  fiducials  were  tracked  by  motion  picture  cameras.  Triangu¬ 
lation  gave  the  sled  fiducials  as  a  function  of  time  in  inches.  The 
acceleration  and  velocity  of  the  sled  relative  to  the  track  were  measured 
in  g  and  ft/s,  respectively. 

As  the  sled  accelerated  down  the  track  away  from  its  rest  position 
under  the  impetus  provided  by  a  piston,  the  subject's  constrained 
yielding  torso  was  dragged  along,  as  was  the  subject's  head  attached 
to  the  torso  by  the  neck.  The  problem  is  to  derive  the  motion  of  the 
head  (translational  and  angular  position,  velocity,  and  acceration) 
by  processing  all  the  data  simultaneously  through  the  Kalman  filter. 


The  nominal  error  standard  deviations  assumed  for  the  experimental 
measurements  are  listed  in  Table  8-1.  The  observation  covariance  matrix 


R  is  diagonal  with  diagonal  elements  the  squares  of  the  standard 
deviations.  Each  observable  has  its  own  units.  The  theoretical 
values  of  the  observables  in  the  units  internal  to  the  software  (cm 
and  s)  are  converted  to  the  observable  units  by  multiplying  by  the 
appropr iate  factors. 


Table  8-1.  Measurement  standard  deviations. 


OBSERVABLE  TYPE 

STANDARD  DEVIATION 

. 

Linear  Piezoresistive 

Accelerometer 

0.1  g 

Photographic 

0.1  inch 

Track  velocity 

0.1  ft/s 

Track  Acceleration 

0.1  g 

The  interval  between  accelerometer  measurements  is  one  millisecond, 
whereas  photographic  measurements  are  two  milliseconds  apart.  An  impact 
event  lasts  about  0.3  seconds.  The  Kalman  filter  allows  different  numbers 
of  observables  at  each  one  millisecond  time  point,  with  the  state  being 
propagated  between  time  points  by  the  dynamical  equations.  The  dynamical 
noise  Q  chosen  in  Section  7.4  for  these  experiments  allows  the  proper 
weighting  of  present  observable  information  relative  to  the  past 
observable  and  initial  condition  information  contained  in  the  propagated 
state  variables. 

8.1  LATERAL  IMPACT  TESTS 

A  diagram  of  the  9-TAP  is  given  in  Figure  8-1.  The  right-hand 
coordinate  convention  of  the  software  is  employed  rather  than  the  left- 
hand  convention  of  the  experiment. 
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There  are  three  accelerometers  near  the  origin  0  of  the  9-TAP 
with  input  axes  along  the  three  coordinate  axes  as  shown.  At  a  distance 
d  =  1.5  inches  =  3.81  cm  along  each  of  the  three  coordinate  axes 

are  two  accelerometers  with  input  axes  perpendicular  to  the  given  axis. 

The  seismic  mass  coordinates  and  input  axis  directions  of  the 
accelerometers  relative  to  the  accelerometer  array  frame  are  given  in 
Table  8-2.  They  are  offset  a  distance  d'  =  0.34  inch  =  0.6  cm 

from  the  center  lines  of  the  array  arms. 

Table  8-2.  Nine  accelerometer  array  coordinates. 

d  =  length  of  array  arms 
=  1.5  inches  =  3.81  cm 

d'  =  offset  of  seismic  masses  from  center  line  of  arms 
=  0.34  inch  =  0.86  cm 


ACCELEROMETER 

POSITION  IN 

ARRAY 

INPUT  AXIS  DIRECTION 

COSINES 

INITIAL  OUTPUT 

WITHOUT  BIAS 

Number 

Name 

*>1 

b2 

*>3 

Cl 

c2 

c3 

B 

x0 

d' 

0 

0 

■fl 

0 

0 

0 

Bb 

*0 

0 

d* 

0 

0 

0 

H 

z0 

0 

0 

-d' 

-1 

-1  g 

B 

Xl 

d’ 

-d* 

0 

1  , 

0 

0 

5 

Zl 

0 

-d 

-d' 

0 

-1 

-1  g 

6 

y2 

-d 

d* 

0 

0 

1 

0 

0 

1 

z2 

-d 

-d* 

0 

0 

-1 

-1  g 

BB 

X3 

d' 

d 

1 

0 

0 

0 

Bfl 

*3 

0 

d 

0 

1 

0 

0 

At  the  initial  time  tg  =  0,  the  sled  is  at  rest  with  zero  velocities 
and  accelerations.  The  track,  sled,  and  accelerometer  array  coordinate 
axes  are  all  initially  parallel  with  zero  relative  velocity  and  accelera¬ 
tion.  The  origin  of  the  sled  frame  coordinates  and  the  position  of  the 


accelermeter  array  frame  relative  to  it  vary  with  each  experiment. 
The  origins  of  the  sled  and  track  frames  are  assumed  to  coincide  at 
the  initial  time. 


8.1.1  Dummy  Test  1363 

Lateral  impact  dummy  sled  test  1363  is  depicted  in  Figure  8-2. 
The  sled  frame  has  origin  in  the  seat,  and  the  initial  position  of  the 
center  of  the  nine-accelerometer  array  and  the  photographic  fiducials 
relative  to  the  sled  frame  are  given  in  Table  8-3. 


Figure  8-2.  Lateral  impact  dummy  sled  test  1363. 
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Table  8-3.  Dummy  lateral  impact  test  1363  initial  accelerometer  array 
and  photographic  fiducial  coordinates  relative  to  the  sled 
frame* . 


9-TAP 

UlO'  z20'  z3o) 


Right  Eye  Fiducial 
l£011'  f012»  foi3) 


Nose  Fiducial 
(f041'  f042'  f043) 


INCHES 

CENTIMETERS 

_ 

_ 

_ 

X 

y 

Z 

X 

Y 

Z 

8.50 

-2.15 

34.80 

21.59 

-5.46 

88.39 

9.00 

-0.65 

32.55 

22.86 

-1.65 

82.65 

10.05 

1.60 

32.15 

25.53 

4.06 

81.66 

Typical  accelerometer  and  photographic  observables  and  Kalman 
filter  post-fit  residuals  are  given  in  Figures  8-3  and  8-4.  The 
residuals  result  from  the  Kalman  filter  trying  to  best  fit  the  data 
from  the  nine-accelerometer  array,  the  sled  accelerometer  and  velocity 
tachometer,  and  the  three  coordinates  of  each  of  the  two  photographic 
fiducials. 

The  Kalman  filter  is  threading  the  states  through  the  observables, 
but  the  residuals  axe  not  as  small  as  might  be  desired.  For  accelerometer 
readings  of  over  10  g,  the  residuals  were  between  -2.5  and  +1.8  g.  The 
photographic  fiducial  residuals  were  between  -1.0  and  +1.2  inches.  The 
sled  acceleration  and  velocity  relative  to  the  track  had  0.2  g  and 
0.5  ft/s  maximum  magnitude  residuals. 

The  dummy  had  a  three-accelerometer  array  at  the  center  of  its 
head.  The  residuals  were  as  large  as  22  g  in  magnitude.  The  most 
probable  cause  of  these  discrepancies  is  some  error  in  the  input 
coordinates  to  the  software. 


Right-handed  coordinate  frame,  so  that  X  is  the  negative  of  the  value 
in  the  left-hand  photographic  measurement  frame. 


z2  OBSERVABLE  (g) 


•»> 


Figure  8-3.  Dunmy  test  1363  accelerometer  observable  and  residual. 
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8.1.2  Human  Test  1313 

Lateral  impact  human  sled  test  1313  is  similar  to  the  dummy  test 
discussed  in  Section  8.1.1.  The  initial  position  of  the  center  of  the 
nine-accelerometer  array  and  the  photographic  fiducials  relative  to  the 
sled  frame  are  given  in  Table  8-4.  The  Kalman  filter  was  again  obviously 
tracking  the  observables,  but  one  of  the  accelerometer  residuals  got 
as  large  as  -3.3  g  and  another  was  as  large  as  +2.8  g  for  accelerometer 
measurements  of  the  order  of  10  g. 


Table  8-4.  Human  lateral  impact  test  1313  initial  accelerometer  array 
and  photographic  fiducial  coordinates  relative  to  the  sled 
frame*. 


INCHES 

CENTIMETERS 

X 

Y 

Z 

X 

Y 

Z 

9-TAP 

(z1Q,  Z20'  z3()) 

6.42 

-1.94 

32.24 

16.31 

-4.93 

81.89 

Right  Eye  Fiducial 
<f011*  f012'  f 0 1 3) 

7.54 

-0.31 

31.92 

19.15 

-0.79 

81.08 

Nose  Fiducial 
(f041'  f042»  f043> 

8.80 

2.08 

32.75 

22.35 

5.28 

83.19 

8.2  PORE-AFT  (EYES  OUT)  IMPACT  TESTS 

Data  from  fore-aft  (eyes  out)  impact  sled  tests  1133  (human)  and 
1229  (dummy)  were  analyzed.  There  was  a  nine-accelerometer  array 
strapped  to  the  subject's  forehead.  In  order  to  simulate  the  situation 
of  a  three-accelerometer  array  afixed  to  the  teeth  inside  the  mouth, 
only  three  of  the  nine  outputs  were  obtained  for  use  in  the  analysis. 


Right-handed  coordinate  frames,  so  that  X  is  the  negative  of  the  value 
in  the  left-hand  photographic  measurement  frame. 
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There  was  photographic  fiducial  data  from  two  head  fiducials 
in  each  experiment.  Only  the  X  and  z  coordinates  versus  time  were 
recorded,  since  the  motion  of  the  subject  was  mainly  in  the  (X,  Z) 
impact  plane.  Therefore,  in  the  Kalman  filter  analysis  the  Y  accelero¬ 
meter  output  was  ignored,  and  only  the  X  and  Z  accelerometer  data 
used  along  with  the  X  and  Z  photographic  fiducial  data. 

In  both  experiments,  the  Kalman  filter  tracked  the  observables 
until  the  start  of  the  impact  event,  at  which  point  the  Kalman  filter 
diverged  from  the  observables,  with  some  residuals  being  greater  than 
103  g.  The  cause  of  this  divergence  requires  further  investigation. 

8.3  ANGULAR  ACCELEROMETER  COVARIANCE  ANALYSIS 

The  Kalman  filter  software  has  the  option  of  reading  a  trajectory 
output  tape  from  a  previous  Kalman  filter  fit,  adding  dummy  observations, 
and  calculating  the  state  covariance  of  the  input  trajectory  for  the 
expanded  observation  set.  The  Kalman  filter  propagation  of  the 
trajectory  is  skipped  and  the  input  trajectory  states  are  unchanged 
with  only  the  state  covariance  being  updated. 

The  output  magnetic  tape  from  lateral  impact  experiment  1363  was 
read  and  a  covariance  analysis  performed  with  the  addition  of  ideal 
angular  accelerometer  observables  in  three  orthogonal  directions.  The 
angular  accelerometer  measurement  error  was  assumed  to  be  0.01  rad/s2. 

The  quaternion  acceleration  covariance  improved  by  almost  an  order  of 
magnitude  versus  not  having  the  angular  accelerometer  observables.  The 
quaternion  velocity  and  position  covariance  improved  by  20  percent  and 
less  than  10  percent  respectively. 

The  covariance  analysis  assumed  that  the  measurements  are  with¬ 
out  bias,  including  dimensional  errors  and  cross  axis  sensitivities  in 
the  accelerometer  array.  The  real  advantage  of  an  angular  accelerometer 
would  be  to  counteract  the  effect  of  such  biases  in  the  piezoresistive 
accelerometer  array.  To  obtain  the  true  effect  of  the  new  observable 
on  the  state  covariance,  bias  states  would  have  to  be  added  to  the 
Kalman  filter  and  their  uncertainty  included  in  the  covariance 
calculations. 
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